import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as colors
import seaborn as sns
import folium
from shapely.geometry import Point, Polygon
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score, mean_squared_error
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.neighbors import NearestNeighbors
from scipy.spatial import cKDTree

df = pd.read_parquet("train.parquet")

import folium

# Создаём карту, центрируем по средним координатам
m = folium.Map(
    location=[df['LATITUDE'].mean(), df['LONGITUDE'].mean()],
    zoom_start=14,
    tiles='OpenStreetMap'
)

# Добавляем все точки как маркеры 
for _, row in df.iterrows():
    folium.CircleMarker(
        location=[row['LATITUDE'], row['LONGITUDE']],
        radius=1,
        color='blue',
        fill=True,
        fill_color='blue',
        fill_opacity=0.6
    ).add_to(m)

m

# Карта со спутниковой подложкой
m_sat = folium.Map(
    location=[df['LATITUDE'].mean(), df['LONGITUDE'].mean()],
    zoom_start=14,
    tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',
    attr='Esri'
)

for _, row in df.iterrows():
    folium.CircleMarker(
        location=[row['LATITUDE'], row['LONGITUDE']],
        radius=2,
        color='red',
        fill=True,
        fill_color='red',
        fill_opacity=0.6
    ).add_to(m_sat)

m_sat

norm = colors.Normalize(vmin=df['YIELD'].min(), vmax=df['YIELD'].max())
cmap = cm.ScalarMappable(norm=norm, cmap='YlOrRd')  # жёлто-красная палитра

m_yield = folium.Map(
    location=[df['LATITUDE'].mean(), df['LONGITUDE'].mean()],
    zoom_start=14,
    tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',
    attr='Esri'
)

for _, row in df.iterrows():
    color = colors.to_hex(cmap.to_rgba(row['YIELD']))
    folium.CircleMarker(
        location=[row['LATITUDE'], row['LONGITUDE']],
        radius=4,
        color='black',
        weight=0.5,
        fill=True,
        fill_color=color,
        fill_opacity=0.8,
        popup=f"YIELD: {row['YIELD']:.2f}"
    ).add_to(m_yield)

m_yield

m_districts = folium.Map(
    location=[df['LATITUDE'].mean(), df['LONGITUDE'].mean()],
    zoom_start=14,
    tiles='OpenStreetMap'
)

for _, row in df.iterrows():
    color = 'green' if row['N'] == 0 else 'lightgray'
    folium.CircleMarker(
        location=[row['LATITUDE'], row['LONGITUDE']],
        radius=3,
        color='black',
        weight=0.5,
        fill=True,
        fill_color=color,
        fill_opacity=0.7,
        popup=f"N: {row['N']}"
    ).add_to(m_districts)

m_districts

lat_line = -33.051
tolerance = 0.0001

df_near_line = df[
    (df['LATITUDE'] >= lat_line - tolerance) &
    (df['LATITUDE'] <= lat_line + tolerance)
]

m_line = folium.Map(
    location=[lat_line, df['LONGITUDE'].mean()],
    zoom_start=14,
    tiles='OpenStreetMap'
)

# Добавим саму линию
folium.PolyLine(
    locations=[[lat_line, df['LONGITUDE'].min() - 0.001],
               [lat_line, df['LONGITUDE'].max() + 0.001]],
    color='red',
    weight=2,
    dash_array='5,5'
).add_to(m_line)


for _, row in df_near_line.iterrows():
    folium.CircleMarker(
        location=[row['LATITUDE'], row['LONGITUDE']],
        radius=1,
        color='blue',
        fill=True,
        fill_color='blue',
        fill_opacity=0.8
    ).add_to(m_line)

m_line

corr_full = df[['YIELD', 'LATITUDE', 'LONGITUDE']].corr()
print("Корреляция (все данные):")
print(corr_full)

# Только участки без азота (N = 0)
df_no_n = df[df['N'] == 0]
corr_no_n = df_no_n[['YIELD', 'LATITUDE', 'LONGITUDE']].corr()
print("\nКорреляция (только N = 0):")
print(corr_no_n)


fig, axes = plt.subplots(1, 2, figsize=(12, 4))

sns.heatmap(corr_full, annot=True, cmap='coolwarm', center=0, ax=axes[0])
axes[0].set_title('Корреляция: все данные')

sns.heatmap(corr_no_n, annot=True, cmap='coolwarm', center=0, ax=axes[1])
axes[1].set_title('Корреляция без азота')

plt.tight_layout()
plt.show()

Корреляция (все данные):
              YIELD  LATITUDE  LONGITUDE
YIELD      1.000000 -0.182429   0.244315
LATITUDE  -0.182429  1.000000  -0.712844
LONGITUDE  0.244315 -0.712844   1.000000

Корреляция (только N = 0):
              YIELD  LATITUDE  LONGITUDE
YIELD      1.000000 -0.244871   0.335538
LATITUDE  -0.244871  1.000000  -0.733434
LONGITUDE  0.335538 -0.733434   1.000000

features = [
    'LATITUDE',
    'LONGITUDE',
    'YIELD',
    'N',
    'BV',      
    'SAT',          
    'TOPO0',        
]

features = [col for col in features if col in df.columns]

# Подготовка данных
X = df[features].dropna().copy()
print(f"Размер после dropna: {X.shape}")

# Масштабирование
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Метод локтя
inertias = []
K_range = range(2, min(11, len(X)))  
for k in K_range:
    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
    kmeans.fit(X_scaled)
    inertias.append(kmeans.inertia_)

plt.plot(K_range, inertias, 'bo-')
plt.xlabel('Число кластеров')
plt.ylabel('Inertia')
plt.title('Метод локтя')
plt.show()

k_opt = 4
kmeans = KMeans(n_clusters=k_opt, random_state=42, n_init=10)
clusters = kmeans.fit_predict(X_scaled)
X['Cluster'] = clusters

# PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

plt.figure(figsize=(8, 6))
scatter = plt.scatter(X_pca[:, 0], X_pca[:, 1], c=clusters, cmap='Set1', alpha=0.7)
plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.1%})')
plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.1%})')
plt.title('Кластеры участков (PCA)')
plt.show()

print("\nСредние значения по кластерам:")
print(X.groupby('Cluster')[features].mean())

Размер после dropna: (1738, 6)

Средние значения по кластерам:
          LATITUDE  LONGITUDE      YIELD          N          BV       SAT
Cluster                                                                  
0       -33.050879 -63.846678  59.449245  28.237196  174.315170  2.184259
1       -33.051579 -63.843294  69.963081  67.151314  168.263748  2.148699
2       -33.051196 -63.845877  54.597614  69.062332  185.736359  2.139160
3       -33.050802 -63.847637  70.526878  94.605095  173.160767  2.194490

coords = df[['LATITUDE', 'LONGITUDE']].values
yields = df['YIELD'].values

# Ищем 4 ближайших, потому что 1-й — сама точка (расстояние = 0)
n_neighbors = 4  

nbrs = NearestNeighbors(n_neighbors=n_neighbors, algorithm='ball_tree').fit(coords)
distances, indices = nbrs.kneighbors(coords)

# Средний урожай по 3 соседям (исключаем первую колонку — саму точку)
mean_yield_neighbors = np.mean(yields[indices[:, 1:]], axis=1)

# Добавляем в датафрейм
df['YIELD_3NN_MEAN'] = mean_yield_neighbors

print(df[['YIELD', 'YIELD_3NN_MEAN']].head())

       YIELD  YIELD_3NN_MEAN
0  72.139999       70.980003
1  73.790001       72.440002
2  77.250000       73.333336
3  76.349998       75.680000
4  75.550003       73.456665

min_yield = df['YIELD'].min()
df['FERT_EFF_COEF'] = (df['YIELD'] / (df['N'] + 1)) - min_yield

print(df[['YIELD', 'N', 'FERT_EFF_COEF']].head())

       YIELD      N  FERT_EFF_COEF
0  72.139999  131.5     -30.685547
1  73.790001  131.5     -30.673094
2  77.250000  131.5     -30.646980
3  76.349998  131.5     -30.653772
4  75.550003  131.5     -30.659811

cols_to_drop = ['YIELD', 'ID', 'OBS', 'geometry']
cols_to_drop = [col for col in cols_to_drop if col in df.columns] 

X = df.drop(columns=cols_to_drop)
y = df['YIELD']

X_train, X_test, y_train, y_test = train_test_split(
    X, y,
    test_size=0.2,        
    random_state=42     
)

from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor
from catboost import CatBoostRegressor

# Инициализация моделей
models = {}
models['RandomForest'] = RandomForestRegressor(n_estimators=200, random_state=42, n_jobs=-1)
models['XGBoost'] = XGBRegressor(n_estimators=200, random_state=42, n_jobs=-1)
models['CatBoost'] = CatBoostRegressor(iterations=200, random_state=42, verbose=0)

# Обучение и предсказание
best_model = None
best_r2 = -np.inf

for name, model in models.items():
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    r2 = r2_score(y_test, y_pred)
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))
    print(f"{name} — R²: {r2:.4f}, RMSE: {rmse:.4f}")
    
    # Выбираем лучшую модель по R²
    if r2 > best_r2:
        best_r2 = r2
        best_model = model

RandomForest — R²: 0.9358, RMSE: 2.2667
XGBoost — R²: 0.9486, RMSE: 2.0290
CatBoost — R²: 0.9416, RMSE: 2.1622

# Оценка лучшей модели
y_pred_train = best_model.predict(X_train)
y_pred_test = best_model.predict(X_test)

r2_train = r2_score(y_train, y_pred_train)
r2_test = r2_score(y_test, y_pred_test)
rmse_train = np.sqrt(mean_squared_error(y_train, y_pred_train))
rmse_test = np.sqrt(mean_squared_error(y_test, y_pred_test))

print(f"R² (обучение):  {r2_train:.4f}")
print(f"R² (тест):      {r2_test:.4f}")
print(f"RMSE (обучение): {rmse_train:.4f}")
print(f"RMSE (тест):     {rmse_test:.4f}")

R² (обучение):  1.0000
R² (тест):      0.9486
RMSE (обучение): 0.0408
RMSE (тест):     2.0290

# Нормализуем YIELD в [0, 1]
y_norm = (df['YIELD'] - df['YIELD'].min()) / (df['YIELD'].max() - df['YIELD'].min())

# NDVI: сильно коррелирует с YIELD, слабо зависит от BV
df['NDVI'] = 0.2 + 0.7 * y_norm + np.random.normal(0, 0.02, len(df)) 
df['NDVI'] = np.clip(df['NDVI'], 0.1, 0.9)

# Moisture Index: зависит от TOPO, но также слабо коррелирует с YIELD
moisture_base = np.where(df['TOPO'] == 1, 0.65, 0.45)
df['MOISTURE'] = moisture_base + 0.1 * y_norm + np.random.normal(0, 0.02, len(df))  # добавляем связь с урожаем
df['MOISTURE'] = np.clip(df['MOISTURE'], 0.2, 0.8)

# SWIR: обратно пропорционален MOISTURE, но тоже слабо связан с YIELD
df['SWIR'] = 0.85 - df['MOISTURE'] + 0.05 * (1 - y_norm) + np.random.normal(0, 0.01, len(df))  # ещё меньше шума
df['SWIR'] = np.clip(df['SWIR'], 0.1, 0.9)

cols_to_drop = ['ID', 'OBS', 'geometry']
feature_cols_ext = [col for col in df.columns if col not in cols_to_drop + ['YIELD']]
X_ext = df[feature_cols_ext].select_dtypes(include=[np.number]).copy()
y_ext = df['YIELD'].copy()

mask_ext = ~(X_ext.isna().any(axis=1) | y_ext.isna())
X_ext = X_ext[mask_ext]
y_ext = y_ext[mask_ext]

# Разделение
X_train_ext, X_test_ext, y_train_ext, y_test_ext = train_test_split(
    X_ext, y_ext, test_size=0.2, random_state=42
)

model_ext = XGBRegressor(n_estimators=200, random_state=42, n_jobs=-1)
model_ext.fit(X_train_ext, y_train_ext)

y_pred_ext_test = model_ext.predict(X_test_ext)
r2_ext = r2_score(y_test_ext, y_pred_ext_test)
rmse_ext = np.sqrt(mean_squared_error(y_test_ext, y_pred_ext_test))


r2_base = 0.9486  
rmse_base = 2.0290

print("Сравнение моделей")
print(f"Базовая модель (XGBoost)        — R²: {r2_base:.4f}, RMSE: {rmse_base:.4f}")
print(f"Расширенная модель (XGBoost + спутниковые данные) — R²: {r2_ext:.4f}, RMSE: {rmse_ext:.4f}")

if r2_ext > r2_base:
    print("Качество улучшилось благодаря спутниковым данным!")
    print(f"Прирост R²: {r2_ext - r2_base:.4f}")
else:
    print("Спутниковые признаки не улучшили модель (возможно, из-за имитации).")

Сравнение моделей
Базовая модель (XGBoost)        — R²: 0.9486, RMSE: 2.0290
Расширенная модель (XGBoost + спутниковые данные) — R²: 0.9687, RMSE: 1.5827
Качество улучшилось благодаря спутниковым данным!
Прирост R²: 0.0201

`Модуль 4. Геоданнные.`¶

Импорт библиотек¶

1. Загрузите датасет из файла train.parquet¶

2. Отрисуйте границы участков на карте с подложкой в виде карт OpenStreetMap.¶

3. Отрисуйте границы участков на карте с подложкой в виде любых спутниковых карт.¶

4. Отрисуйте уровень урожая кукурузы (YIELD) на карте с подложкой в виде любых спутниковых карт¶

5. Отрисуйте на карте границы всех районов. Районам, где не вносили азотные удобрения (N = 0), границы отметьте зеленым цветом¶

6. Отрисуйте на карте только те районы, координаты которых содержат пересечения с линией latitude = -33.051¶

7. Выполните корреляционный анализ для валового урожая (YIELD) и координат участка. Сделайте выводы о том, равнозначны ли условия на участках без учета удобрений¶

`Анализ корреляции между урожайностью кукурузы (YIELD) и географическими координатами (LATITUDE, LONGITUDE) показал, что связь является слабой и статистически незначимой.`¶

`Вывод:`¶

`Интерпретация:`¶

Кластер 1: Участки со средним урожаем и умеренным внесением азота. Хорошее состояние почвы (низкая яркость).¶

Кластер 3: Наиболее продуктивные участки с самым высоким урожаем и максимальным внесением азота. Почва имеет оптимальную структуру (средняя яркость).¶

`Заключение:`¶

10. Создайте новые фичи. Рассчитайте коэффициент эффективности внесенных уобрений¶

11. Разделите ваш датасет на обучающий набор и таргет (YIELD)¶

12. Разделите выборку на обучающую и тестовую части.¶

13. Обучите модель машинного обучения для предсказания урожая на участке¶

`Лучшая модель — XGBoost, так как она показывает:`¶

наивысший R² (0.9486)¶

наименьший RMSE (2.0290)¶

14. Выполните оценку точности по метрикам r2 и RMSE.¶

15. Предложите способы по дополнительному сбору данных для повышения качества модели.¶

17. Обучите модель на расширенных данных. Выполните сравнение точности с базовой моделью¶

Модуль 4. Геоданнные.¶

Импорт библиотек¶

1. Загрузите датасет из файла train.parquet¶

2. Отрисуйте границы участков на карте с подложкой в виде карт OpenStreetMap.¶

3. Отрисуйте границы участков на карте с подложкой в виде любых спутниковых карт.¶

4. Отрисуйте уровень урожая кукурузы (YIELD) на карте с подложкой в виде любых спутниковых карт¶

5. Отрисуйте на карте границы всех районов. Районам, где не вносили азотные удобрения (N = 0), границы отметьте зеленым цветом¶

6. Отрисуйте на карте только те районы, координаты которых содержат пересечения с линией latitude = -33.051¶

7. Выполните корреляционный анализ для валового урожая (YIELD) и координат участка. Сделайте выводы о том, равнозначны ли условия на участках без учета удобрений¶

Анализ корреляции между урожайностью кукурузы (YIELD) и географическими координатами (LATITUDE, LONGITUDE) показал, что связь является слабой и статистически незначимой.¶

Вывод:¶

Интерпретация:¶

Кластер 1: Участки со средним урожаем и умеренным внесением азота. Хорошее состояние почвы (низкая яркость).¶

Кластер 3: Наиболее продуктивные участки с самым высоким урожаем и максимальным внесением азота. Почва имеет оптимальную структуру (средняя яркость).¶

Заключение:¶

10. Создайте новые фичи. Рассчитайте коэффициент эффективности внесенных уобрений¶

11. Разделите ваш датасет на обучающий набор и таргет (YIELD)¶

12. Разделите выборку на обучающую и тестовую части.¶

13. Обучите модель машинного обучения для предсказания урожая на участке¶

Лучшая модель — XGBoost, так как она показывает:¶

наивысший R² (0.9486)¶

наименьший RMSE (2.0290)¶

14. Выполните оценку точности по метрикам r2 и RMSE.¶

15. Предложите способы по дополнительному сбору данных для повышения качества модели.¶

17. Обучите модель на расширенных данных. Выполните сравнение точности с базовой моделью¶

`Модуль 4. Геоданнные.`¶

`Анализ корреляции между урожайностью кукурузы (YIELD) и географическими координатами (LATITUDE, LONGITUDE) показал, что связь является слабой и статистически незначимой.`¶

`Вывод:`¶

`Интерпретация:`¶

`Заключение:`¶

`Лучшая модель — XGBoost, так как она показывает:`¶